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Abstract 

A toy model of two dimensional nanoindentation in finite crystals is proposed. 
The crystal is described by periodized discrete elasticity whereas the inden- 
ter is a rigid strain field of triangular shape representing a hard knife-like 
indenter. Analysis of the model shows that there are a number of discontinu- 
ities in the load vs penetration depth plot which correspond to the creation 
of dislocation loops. The stress vs depth bifurcation diagram of the model 
reveals multistable stationary solutions that appear as the dislocation-free 
branch of solutions develops turning points for increasing stress. Dynamical 
simulations show that an increment of the applied load leads to nucleation 
of dislocation loops below the nanoindenter tip. Such dislocations travel 
inside the bulk of the crystal and accommodate at a certain depth in the 
sample. In agreement with experiments, hysteresis is observed if the stress 
is decreased after the first dislocation loop is created. Critical stress values 
for loop creation and their final location at equilibrium are calculated. 
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1. Introduction 



Regularizations of continuum theories typically resolve the difficulties of 
the latter at small scales and often allow precise calculations and descrip- 
tions at these scales. This is the case of the Navier-Stokes equations which 
regularize the inviscid Euler equations of Fluid Mechanics (ll|), of the lat- 
tice regularizations of Quantum Field Theory and phase transitions ([2|), etc. 
Classical elasticity theory is unable to describe the nucleation and motion 
of crystal defects (^, dislocations (4), cracks (5j 6), phase boundaries 



111 ). Lattice regu- 



or more complex phenomena such as friction ((8|; |9|; lid : 
larizations of elasticity can describe the structure and motion of nucleated 
defects and can be analyzed to extract qualitative and quantitative informa- 



tion about t 



standing (JIO 



ie_ phenomena at hand so as to provide a deep physical under- 
. In the case of dislocations, periodized discrete elasticity can 



describe depinning of dislocations at the Peierls stress (Il2l ). dislocation cores 



and dislocation interaction (13), stable defects corresponding to dis^ 



in graphene membranes and instability of Stone Wales de: 



homogeneous nucleation of dipoles in a sheared lattice (Il6l ) 



'ects (jjj; 



ocations 



ISl). and 



Metals usually contain a great number of dislocations whose motion, cre- 
ation, annihilation and interaction are largely responsible for plastic behav- 
ior (jj). Introducing defects in a crystal typically impedes dislocation motion 
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and multiplication thereby strengthening the material (strain hardening). 
At small length scales (such as those intervening in compression of thin 
whiskers), dislocations may leave the sample which results in its harden- 
ing via dislocation starvation (jlTI ) as observed in compression of nanopillars 



(jlSl ). Incipient plasticity occurs when defects are created in a hitherto perfect 



crystal by different means. Nanoindentation experiments are excellent ways 



to probe incipient plasticity (19 



2C 
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22I ) and so are indentation experi- 



ments in colloidal crystals (1231 ). In these experiments, the penetration depth 
inside the crystal is measured as the load on the indenter increases, which 
results in a discontinuous load versus depth diagram. The discontinuities in 
the diagram are thought to indicate dislocation nucleation inside the crystal. 
Different types of calculations have been used to interpret nanoindentation, 



ranging from atomistic simu 
or combinations thereof 



ations to continuum mechanics interpretations 



21 



26|). 



In this paper we present and analyze an atomistic toy model of nanoinden- 
tation to show that discontinuities in the load vs penetration length diagram 
do correspond to nucleation of dislocation loops q We also use this model 
and its analysis to find different loading-unloading stress vs depth curves, in 
agreement with experimental observations (fiol ). We have used the AUTO 
software (1271 ) to find and numerically continue the branches of stationary 



^In a 2D model, dislocation dipoles are pairs of edge dislocations of opposite Burgers 
vectors with infinitely long extra half rows of atoms that go from the dislocation point 
outwards. Dislocation loops are pairs of edge dislocations with opposite Burgers vectors 
sharing a finite segment of extra atoms that join their respective dislocation points. See 
Figures 10 and 9 of (jjj). 
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configurations of our model starting from the perfect lattice. Thus we are 
able to calculate both stable and unstable stationary configurations and the 
bifurcation diagram (maximum strain/dimensionless stress F vs penetration 
depth 6), which is something that cannot be done using Molecular Dynam- 
ics (MD) simulations due to neighbor upgrading protocols. The bifurcation 
diagram exhibits multistability of stationary solutions whose configurations 
display an integer number of dislocation loops. Starting from the unstressed 
crystal and as F increases, S increases until a turning poinj^ is reached. As 
F surpasses a critical value Fc, a dislocation loop is nucleated and the profile 
evolves to that of the next stable solution branch in the bifurcation diagram 
(having a larger 6). This gives rise to a discontinuity in the load vs depth 
diagram for static solutions. Further increments of the load give rise to new 
discontinuities when successive turning points of the different stable solution 
branches are reached. In this simple model, nucleation of dislocation loops 
is a first order phase transition and hysteresis is possible: for F > Fc, there 
is an abrupt jump of 6 due to dislocation loop nucleation and the unloading 
curve follows the second stable branch which has larger 6 for the same F as 
in the loading curve. For our simple model, the second stable branch ends at 
a positive value of F. Further decrease thereof provokes reabsorption of the 
nucleated dislocation loop at the indenter tip and a sudden decrease of 6 to 
its value on the elastic curve. 

Our toy model is as follows. Consider a 2D simple cubic lattice with 
displacement vector measured in units of the lattice constant I, lattice points 



"Points with zero or infinite slope where a branch of solutions changes from stable to 



unstable are turning points. See Fig. II. 1 of (1281 ) 
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Figure 1: Modeling of a nanoindcntation experiment. The stress field scales with the 
numerical continuation parameter F, such that the maximum downward vertical displace- 
ment at the central part of the upper surface satisfies D2 viq^/2^Ny — ■ 



labelled by indices (i, j), i = 1, . . . , iVj, and j = 1, . . . , Ny. The nanoinden- 
ter at the central region of the upper surface is a hard 'knife' represented 
by a fixed strain field which decreases linearly from zero to —F and then 
increases to zero again, as depicted in Fig. [H At the other boundaries, the 
displacement vector is zero. The indenter may cause gliding of atom columns 
in the vertical direction and therefore the displacement vector changes more 
in the vertical direction than in the horizontal one. Then a simple dis- 
placement vector {0,Vij) captures the qualitative features of the physical 
system. It is convenient to use coordinates {x,y) with x = i — {N^ + l)/2 
and y = j — {Ny + l)/2 whose origin is the lattice center. The displacements 
obey the following nondimensional equations: 

^~d^ + = ~ + Vij-i + A [ga{vi+ij - Vij) + ga{vi-i,j - Vij)], (1) 

where ga{x) is a one-parameter family of periodic functions of x with period 
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2a J sin (f ) , -a<x< a, 

^ I sm ( ^ 1 , a < X < 1 — a, 

and < a < 1/2, such that ga{x) ~ x as x — > (see below for motivation). 
The boundary conditions are 

D2Vt,Ny = Vi,Ny - Vi^Ny-1 = /(^^ 

= 0, if either j = 1, i = 1, or i = N^- (4) 

For a symmetric and centered sharp indenter as depicted in Fig. [1] with an 
even number of atoms, w, and for even N^-, we have 



w 



i 2 ' 2 — ^ — 2 ' 

^ + 1-^, f + i<z<te, (5) 



0, otherwise. 

This indenter has a rectangular cross section S = Lwl, vol <^ L, where / is the 
lattice constant and L is the length of the knife. A more realistic modeling 
of the knife would require considering horizontal displacements which we are 
ignoring in our model. Note that f{Nx/2) = /(I + Nx/2) = 1 and that 
J2ifi'^) = 1 +w/2. In Eq. ([T]), A = C44/C11 provided we consider cubic 
crystals with elastic constants Cu, Cu, C44. The vertical component of the 
stress tensor is simply Cu -D^f in our model, and therefore the strain at 
the surface given by Eq. ([3]) is also the nondimensional applied stress cr22/Cii. 
We also have cr2i/Cii = A ga{D^Vi j)/2. The relation between the stress at 
the surface and an apphed load P is 

P = — ^sinLp{i) a2iii, Ny) ^ cos 0-22 (i, iVj/), 

^ i=2 ^ i=2 



-8 -6 -4 -2 2 4 6 
x/l 



Figure 2: An edge dislocation parallel to the z axis with Burgers vector (0,-1). 

where ip{i) is the angle between the normal to the indenter surface and the 
y axis, with tany^ = D^Vi ^y- We obtain 

A = TT^ = — y'cos(/?(i)/(i) + — Vsin(/?(2)^a(L)i-T;i,jvJ, (6) 

1=2 1=2 

which relates the nondimensional load A to the stress F. If we select a nondi- 
mensional time scale Cut/{pP'y) — > t, then a = 1 and m = Cii/(p/^7^), 
where 7 is a friction coefficient with units of frequency, p is the mass density 
and Vij is measured in units of /. With this choice of scales, we can consider 
the overdamped case with m = 0. On the other hand, if we select a nondi- 
mensional time scale C\{^t/{lp^/'^) — > t, then m = 1 and a = ^7a/ p/Cu. 
With this second choice of scales, we can consider the conservative case with 
a = 0. 

If ga{x) = X in Eq. ([1]), we obtain discretized scalar linear elasticity. 
Why do we have to use the periodic functions ([2])? To allow atoms change 
neighbors during dislocation motion. The fact that the functions ([2]) are 



periodic allows the atoms to change neighbors while the computational grid 
remains unchanged (and so the dynamical system ([1]) whose linear stability 
will be analyzed). 

Let us explain this idea. Assume that we have an extra half row of atoms 
in the square lattice that is parallel to the positive x direction (edge disloca- 
tion with Burgers vector (0, —1)), as depicted in Fig. [2l If these atoms move 
one step downwards and occupy the equilibrium positions of their neighbor- 
ing atoms, the latter are displaced downwards and become the extra half row 
representing the edge dislocation that has therefore moved one step down- 
wards. If these new atoms move one step downwards in the same way, a new 
half row replaces them and the dislocation has moved another step. Dur- 
ing this glide motion, the leftmost atom in the half row changes neighbors 
horizontally and only once whereas all the other atoms in the half row con- 
tinue having the same neighbors. Even though the dislocation can glide for 
a long distance, single atoms move very little and only a very small fraction 
of these atoms change neighbors. Instead of updating neighbors during glide 
motion, we can construct a periodized discrete elasticity model by using a 
nonlinear periodic function in ([1]), with period equal to the lattice space 
and ga{x) ~ x as a; — ^ 0. A simple example is ([2]). This model changes 
neighbors without updating, thereby allowing atomic half rows to glide in 
the y-direction. 

As the applied stress becomes sufficiently large, there appear edge dis- 
location loops whose Burgers vectors are directed along the y axis. Gliding 
along other directions is not possible in this model: we would need a two- 
component displacement vector and a periodic function of discrete differences 
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along the x axis (jl3l ). The parameter a controls the asymmetry Qa- a in- 
creases, the interval over which g'^ix) > increases at the expense of the 
interval over which the slope of Qa is negative. As a increases, the (dimen- 
sionless) Peierls stress ap needed for a dislocation to start moving increases 



whereas the size of the dislocation core and its mobility both decrease ( 113 



The value of a can be selected so that the Peierls stress calculated from (1) 
fits values measured in experiments or calculated using MD. Namely, for a 
16 X 30 lattice with A = 0.2258 (corresponding to gold) we get ap = 0.004 
for a = 0.2 compared to cxp = 0.03 for a = 0.4. 

In the symmetric case a = 1/4, ([1]) and ([2]) are the governing equations 
of the interacting atomic chains model (1291 ). 



2. Methodology and Results 

We consider parameters A = 0.2258 (gold), a = 0.2, for a 22 x 42 lattice 
and a four atom indenter, w = 4. Similar results were obtained when the 
same indenter is placed in the center of the upper surface of other lattices 
with even Nx and Ny. Asymptotically stable stationary solutions of the model 
equations with boundary conditions ([3]) - (jl]) are obtained by numerically 
solving the equations with appropriate initial conditions and overdamped 
dynamics, m = 0, a = 1. For conservative dynamics, m = 1, a = 0, the 
same stationary solutions are stable but not asymptotically stable. Unstable 
solutions are unstable no matter which dynamics is considered. 

At zero applied load the perfect lattice is found. As the dimensionless 

■^Note that the parameter a used in ([j3) corresponds to —a + 1/2 in ([2]) and therefore 
the Peierls stress in Figure 2 of |13i) decreases as a increases. 
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Figure 3: (a) Maximum vertical nondimensional stress F vs nondimensional indenter 
penetration depth 6 obtained by numerical continuation of the stationary solution at 
F = (5 = for a 16x42 lattice, with A = 0.2258 (gold), a = 0.2. By increasing adiabatically 
F from zero, we sweep this bifurcation diagram following the path A A' B ^ B' ^ 
C ^ C ^ D ^ D' ^ E ^ E' . (b) Similar diagram of nondimensional applied load 
vs nondimensional penetration depth. Only the long stable portions of the bifurcation 
diagram arc shown. 
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applied load F increases from zero, we can monitor the penetration depth of 
the indenter, 

^ iN^+w)/2 
i=l+(N^-w)/2 

and obtain the relation between F and 5. We have used the AUTO software 



( 1271 ) to find and numerically continue the branches of stationary solutions of 



our model starting from the perfect lattice with F = 0, 5 = 0. The results 
are represented in Fig. [3k, whereas Fig. [3b shows the corresponding load vs 
penetration depth diagram (in dimensionless units). The bifurcation diagram 
shows a connected branch of stationary solutions which begins at the origin 
(marked with A in Fig. ^i). Each point of the long stable portions of this 
branch corresponds to a single configuration of the lattice. There are many 
turning points (local maxima or minima, some of them indistinguishable at 
the scale of the plot) connecting stable and unstable regions. The main 
turning points point out to the occurrence of different nucleation events, 
which we shall now describe. 

The longest stable portions of the stationary branch end at turning points 
A\ B', C, D' and E' in increasing order of 5 (solid lines in Fig. The cor- 
responding configurations have 0, 1, 2, 3 and 4 dislocation loops, respectively, 
as shown in Figs. [1^, c and e. We may call AA' the dislocation-free elastic 
branch, whereas the other long stable portions of Fig. [3k signal the beginning 
of plasticity due to loop formation. The component ei,2 = ga{vi+ij — Vij)/2 
of the strain field (proportional to the dimensionless shear stress) provides 
a good visualization of the dislocation loops as Figures [Dd, d and f show. 
Each loop comprises two dislocations whose cores are located at points with 
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Figure 4: Stable configurations at the turning points (a) A', (c) B' and (e) C" of the 
bifurcation diagram in Fig. [3] Panels (b), (d) and (f) visualize the component ei.2 = 
ga{vi+ij — Vij)/2 of the strain tensor corresponding to A' , B' and C", respectively. 
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2 = 1 + {Nx — w)/2 and i = {N^ + w)/2 below the indenter. Using the 
lattice constant as unit of length, the first dislocation is at the column 
i = 1 + [Nx — w)/2 and it has Burgers vector (0,-1), whereas the sec- 
ond dislocation is located at the column i = [N^ + w)/2 and it has Burgers 
vector (0, 1). 

Note that after some of its turning points the solution branch bends over 
itself, and connects back to earlier limit points (C, D', E' in Fig. [3^) by 
means of unstable portions. To each point of the short stable portions found 
after these limit points, there correspond two lattice configurations, which 
are symmetrical with respect to x = 0. They contain one extra edge disloca- 
tion in addition to the dislocation loops that the configuration corresponding 
to the preceding long stable portion with smaller 6 may contain. In one con- 
figuration, the dislocation point of the extra dislocation with Burgers vector 
(0, 1) is located at i = {Nx + w)/2, at some distance below the indenter (as 
occurs in the short stable portions after A' and B'). The other configura- 
tion has a dislocation point on the same row, with opposite Burgers vector 
(0, —1), which is located at z = 1 + [N^ — w)/2. As it will be explained later, 
these short stable portions are not reached by any of our dynamical exper- 
iments. However, if the configurations corresponding to unstable branches 
are depicted, nucleation or disappearance of the extra edge dislocation is ob- 
served as we run AUTO along the unstable branch towards the turning point 
where it joins a stable branch. 
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3. Adiabatic sweeping of the bifurcation diagram 

By increasing adiabatically F from F = (up-sweep), the stable portions 
AA\ BB', CC, DD', EE' of the solution branch in the bifurcation diagram 
of Fig. [3^ are successively swept. A' marks the critical value F^ above which 
dislocation loops are nucleated. B', C, D' and E' are the last points found 
in the corresponding stable portion of the solution branch before a new dislo- 
cation loop is nucleated and a jump to another stable portion occurs during 
adiabatic up-sweep. At these points, with stresses F^/, Fb', Fc, F^/ and so 
on, 5 jumps to the corresponding value in the next branch, as indicated by 
the arrows in Fig. and a new dislocation loop is formed. What is ob- 
served during slow up-sweep by numerically solving ([T])? After F surpasses 
Fa' in Fig. [3^ with configuration as in Fig. H^, a dislocation loop is nucleated 
immediately below the indenter tip and it glides downwards until it reaches 
its stable position inside the sample at y/l = —10 (configuration B). The 
jump A' ^ B in Fig. [3^ has taken place. By adiabatically increasing F, we 
may reach B' (cf. Fig. Ht and d), which is almost identical to B. After F sur- 
passes Fb', a new dislocation loop is nucleated at the indenter tip, and now 
both loops glide downwards until they reach their respective stable positions 
y/l = 0, —14 in configuration C. Further increase of F leads to configuration 
C (see Fig. H^). Additional up-sweep repeats this pattern: old dislocation 
loops glide downwards and new ones appear beside the indenter tip. 

If we start relaxing adiabatically the load on the nanoindenter after F 
has reached Fb, then the displacement S does not go back to the value Sa'- 
Instead, its value decreases following the stable branch containing the point 
{6b, Fb) in Fig. ^1^. As F decreases, the loop created in the middle of the 
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sample glides upwards towards the indenter until the turning point with 
lowest value of F is reached. At this stable configuration, the loop is located 
at y// = 5. An additional decrease of the load results in a downward jump to 
the elastic solution branch. During this dynamical process, the dislocation 
loop glides upwards and disappears as it reaches the indenter (in 3D loops 
may be anchored due to cross-slip or to preexisting defects (|3|)). In a standard 
model of plastic behavior, we would expect that the second stable branch of 
stationary solutions extend to F = 0, 5 > without losing the dislocation 
loop. 

There are other stable portions of the stationary solution branch in Fig. [3^ 
that cannot be reached by up- or down-sweeping the bifurcation diagram: the 
short portions with intermediate values of 5 between those of two long stable 
portions. 

4. Results for other parameter values 

For the same values of a and A, we have changed the position of the 
indenter, made it wider, changed the size of the lattice or considered odd 
Nx- In all these cases, the F (or A) vs 5 diagram is similar to that in Fig. [3l 
as shown in Fig. [S] for a lattice of different size. The unstable portions of 
the bifurcation diagram change although not the long stable portions, and 
therefore the responses to adiabatic up- and down-sweeping the diagram are 
the same. For parameter values corresponding to the point 5', the dislocation 
loop is located at a height roughly two thirds of However there are some 



*For = 16, 5B'/Ny = 0.625, 0.661 and 0.668 for Ny = 30, 42 and 58, respectively. 
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differences worth noticing. For a 22 x 42 lattice, increasing the indenter width 
from w = 4to w = 10 increases only slightly Fa', from 0.1575 to 0.1598, but 
now the cavities created by the indenter are not rectangular as in Fig. |H 
Instead, they exhibit steps which made them triangular in shape. Similar 
results are found for odd A^^^ when the maximum f{i) is reached at only 
one atom, i = [Nx/2]. For a 17 x 30 lattice with w = 5, the bifurcation 
diagrams are similar and triangular cavities are observed. A shorter indenter 
with w = A placed on atoms 7, 8 , 9 and 10 of a 17 x 30 lattice again 
gives rise to bifurcation diagram similar to Fig. [5] but now some of the short 
intermediate stable portions having one extra unpaired edge dislocation can 
be distinguished in the diagram because their corresponding configurations 
are no longer symmetric. Actually, it is not always after the third turning 
point (C in Figs. [3] and [5]) that these short intermediate stable portions 
correspond to two solutions. We find that for a 22 x 22 lattice it occurs 
after the second turning point, B'. If we use A = 0.448 (corresponding to 
copper), for a 16 x 30 lattice the results are qualitatively similar to those so far 
presented except that: (i) the point A' has a larger F = 0.2191 and a smaller 
S = 0.7958 (compared to F = 0.1553 and 6 = 0.8128 for gold in Fig. El] and 
(ii) close to the minima of the bifurcation diagram, the branch of stationary 
solutions is wiggly with several turning points separating alternating stable 
and unstable solutions. These wiggly portions are also present in Fig. [5] but 
they are not visible at the scale chosen in the figure. 

As mentioned before, augmenting a in Eq. ([2]) increases Fa' and the 
Peierls stress and reduces the size of defect cores. For a 16 x 30 lattice with 
w = 4, A = 0.2258 (gold) and a = 0.4, we get Fa' = 0.2558, 6a' = 1.2283 
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(compared to F = 0.1553 and S = 0.8128 for a = 0.2) and it is quite hard 
to move dislocations (the dimensionless Peierls stress is Up = 0.03, compared 
to cTp = 0.004 for a = 0.2). The bifurcation diagram changes in that the 
number of turning points increases enormously. However, the stable portions 
of the diagram have the same configurations and interpretations as those in 
Figs. [3] and [51 Then up-sweeping the diagram produces the same sequence 
of discontinuities due to formation of dislocation loops. 

5. Conclusion 

We have presented a toy model of 2D nanoindentation and incipient 
plasticity in a dislocation-free crystal based on periodized discrete elasticity. 
Analysis and numerical simulation of the model confirm that the discontinu- 
ities in the diagram of stress (or load) vs penetration depth are associated 
with the nucleation of dislocation loops (if stress is adiabatically increased) or 
with their absorption in contact with it (if stress is adiabatically decreased). 
Our simulations and analysis of the bifurcation diagram of solutions show 
that hysteresis occurs, leading to unloading curves which are different to the 
loading ones, in agreement with experimental observations. The analysis 
of the results obtained with AUTO yields the critical stress for dislocation 
nucleation, the type thereof and the depth at which dislocation loops are 
accommodated in their equilibrium positions. 
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